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ABSTRACT 

We have modelled the observed color-magnitude diagram (CMD) at one location in 
M33's outskirts under the framework of a simple chemical evolution scenario which 
adopts instantaneous and delayed recycling for the nucleosynthetic products of Type 
II and la supernovae. In this scenario, interstellar gas forms stars at a rate modulated 
by the Kennicutt-Schmidt relation and gas outflow occurs at a rate proportional to the 
star formation rate (SFR). With this approach, we put broad constraints on the role 
of gas flows during this region's evolution and compare its [cc/Fe] vs. [Fe/H] relation 
with that of other Local Group systems. We find that models with gas inflow are 
significantly better than the closed box model at reproducing the observed distribution 
of stars in the CMD. The best models have a majority of gas inflow taking place in 
the last 7 Gyr, and relatively little in the last 3 Gyr. These models predict most stars 
in this region to have [cr/Fe] ratios lower than the bulk of the Milky Way's halo. 
The predictions for the present-day SFR, gas mass, and oxygen abundance compare 
favorably to independent empirical estimates. Our results paint a picture in which 
M33's outer disc formed from the protracted inflow of gas over several Gyr with at 
least half of the total inflow occurring since z ~ 1. 

Key words: galaxies: abundances, galaxies: evolution, galaxies: individual: Messier 
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1 INTRODUCTION 

, The distribution of stars in a color-magnitude diagram 
(CMD) is a sensitive function of the past and present star 
formation rate and chemical composition. This fact allows 
the extraction of the star formation history (SFH) and chem- 
ical enrichment history (CEH) of a galaxy by fitting its ob- 
served CMD with a model CMD whose SFH and CEH are 
known. The results of this synthetic CMD fitting method can 
ultimately yield insights into the processes shaping galaxy 
evolution, such as gas flows, energetic feedback, satellite ac- 
cretion, tidal interactions, and dynamical mixing. Never- 
theless, this method is only one of many employed in the 
domains of near-field cosmology a nd galactic paleontology 
|Freeman fc Bland-Hawthornll20ul '). 

Another, equally important method, is chemical mod- 
elling, whose primary goal is to reproduce the elemental 
abundance distribution in a galaxy's stars. This method in- 
volves tracking the effects of gas flows, star formation, and 
stellar nucleosynthesis on the abundances and abundance ra- 
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tios of certain elements. For example, the abundance of the 
a-elements (O, Ne, Mg, Si, S, Ca, Ti) relative to iron, com- 
monly plotted as [a/Fe] vs. [Fe/H], is an extremely useful 
tool to discriminate between different models for a galaxy's 
evolution because the a-elements and iron have different 
production sites and time-scales. The a-elements are pro- 
duced mainly in the hydrostatic burning phases and explo- 
sive deaths (supernovae Type II, or SNe II) of massive stars 
(M > 8 Mq), which have short lifetimes (< 10 7 yr). The 
majority of iron, on the other hand, is created in supernovae 
Type la (SNe la), which are the thermonuclear explosions of 
carbon-oxygen white dwarfs in binary systems and, there- 
fore, typically occur on a much longer time-scale of ~ 1 Gyr 
(|Greggidl20u5l h 

The methods of CMD fitting and chemical modelling 
are complementary, but they have evolved mostly indepen- 
dently. Some examples in the recent literature which com- 
bine the two methods generally use the results of CMD fit- 
ting as e xternal inputs to che mical evolution models. For 
example, lAparicio et alj (Il997l ) derived the first SFH of the 
Local Group (LG) dwarf galaxy, LGS 3, from CMD fit- 
ting and used the chemical evolution equations to see how 
much g as inflow and outfl ow was required to match their 
results. ICarigi et all l|2002f ) computed the chemical evolu- 
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tion of four dwarf spheroidal (dSph) satellites of the Milky 
Way (MW) . As external constraint s, these authors u s ed the 
SFHs of these galaxies deriv ed bv iHernandez et al.l (|2000T ) 
from a CMD fitting method. iDolphin l (|2002| ) calculated the 
SFHs of 6 dSph MW satellites using synthetic CMD fitting, 
and then lLanfranchi fc Matteuccil l|2003l . l2004h adopted his 
SFHs as inputs to their chemical evolution models for the 
same galaxies. Similarly, iFenner et ail (|2006T l modelled the 
chemica l evolution of the Scu lptor dSph using the SFH de- 
rived bv lDolphin et all (|2005T l from a CMD analysis . Lastly , 
after adopting the SFH derived by IWvderl (|200ll . l2003h . 
ICarigi et af] (|2006r i modelled the evolution of NGC 6822 in 
a cosmological framework. 

Conversely, the results of chemical modelling can be 
used as external inputs t o CMD fitting. For example, 
IPagel fc Tautvaisiene! (|l998l . hereafter PT98) modelled the 
chemical evolution of the Large Magellanic Cloud (LMC) 
and Small Magellanic Cloud (SMC), adopting gas inflow and 
non-selective galactic winds (i.e., the wind efficiency was the 
same for all chemical elements). These authors tuned the 
model parameters to match the observed elemental abun- 
dances of clusters and supergiant stars in these systems. 
Some subsequent studies a dopted the LMC's age-metallicity 
relation (AMR) derived by |PT98l to e xtract its SFH from the 
CMD and luminosity fun ction (e.g., iHoltzman et al.lll999l; 
ISmecker-Hane et al. 1 120021 ). 

All the aforementioned studies combined CMD fitting 
and chemical modelling in a two-step process, in which the 
first step was done independently of the second. However, 
the two steps are inextricably linked because an increase in 
the star formation rate (SFR) speeds up the chemical en- 
richment. Gas flows into or out of the system can change 
this coupling making them important ingredients to include 
in any model. Therefore, one drawback of using the CMD- 
derived SFHs as independent constraints on the chemical 
evolution models is that the SFHs are not necessarily phys- 
ically self-consistent under the action of processes like nu- 
cleosynthesis, galactic winds, and gas accretion, processes 
which are fundamental to the chemical models themselves. 
Similarly, using the results of chemical evolution models as 
independent constraints on CMD fitting is not necessarily 
self-consistent because the models have a particular form 
for the SFH, which is the very thing the CMD fitting is 
supposed to derive. 

The recent works of llkuta fc Arimotol (|2002l ) and 
lYuk fc Led £2007) represent a significant improvement 
over the studies mentioned above because they incorpo- 
rate CMD fitting and chemical modelling simultaneously. 
llkuta fc Arimotol ((2002) computed a few closed box evolu- 
tion models for several dSph satellites of the MW. They 
performed a qualitative comparison of their model CMD 
and [Mg/Fe] vs. [Fe/H] relation with what was observed and 
found a reasonable agreement, but they had to invoke gas 
stripping via ram pressure or tidal shocks to reconcile the 
present day gas fraction of their closed box models (~ 97%) 
with the observe d values (~ 0%). Yuk fc Leel (|2007f ) improve 
upon the work o f llkuta fc Arimotol (|2002F by quantitatively 
fitting a closed box chemical model to the CMD of 1C 1613, 
a relatively isolated LG dwarf irregular galaxy. Their model 
SFH and AMR are in good agreement with previous inde- 
pendent determinations based on the canonical CMD fitting 
method, lending support to both the old and new methods. 



Moreover, their predicted present-day oxygen abundance is 
consistent with the observed value. 

There are many lines of evidence that suggest 
all galaxies ev olve as closed boxes. As originally 
pothesized by Larson! (|l974T ) and exemplified by 



llkuta fc Arimotol (2002) results, gas outflows could 
plain the lack of ga s in dSphs despite their 
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metallicities (see also lLanfranchi fc Matteuccil l|2004) ). 
Other evidence for gas outflow s includes the abun- 
dances of metals in the IGM (|Edge fc Stewartl Il99ll : 
iMushotzkv fc Loewensteinl fl997T ). the mass-metallicity re- 
lation of gal axies at low and high redshift (e.g 



Garnettl |2002|; iTremonti et all |2004| ; IPilvugin et~all |2004| ; 
Erb et al.ll2006l). extended extra-planar HI gas in spirals 



( Fraternali et al ] 12004 lOosterloo et al.H2007r ). high velocity 
clouds (HVCs) around the MW with near so l ar metallicity 
(Ivan Woerden fc Wakkenl2004 IWakkedl200ll ; iRichter et all 
l200lT ) , extra-pla nar optical and X-ray emission around star- 
burst galaxies dHeckman et al. 1990l; Lehnert fc Heckmanl 
ll996l ; lMartin et al.ll2002l ; IStrickland et al.ll2004 ), and veloc^ 
ity shifts o f high-ion absor ption lines in damped Lyman- 
q systems llFox et al.ll2007t) and L y man break galaxies a t 
z ~ 3 - 4 (e.g.. |Pettini et al.ll200ll ; Udelberger et al.ll2003r i. 
For a more detailed r eview, we refer the interested reader to 
IVeilleux et all l|2005l ). 

Evidence for gas inflows includes the so-called G- 
dwarf problem, which is the fact that the metallicity 
distribution function of low mass, long lived stars ob- 
served in the Solar vicinity (SV) and in many other 
galaxies is too narrow and contains too few met al poor 
stars c ompared to the closed box m odel (e.g., iTinslevI 
19751; iRocha-Pinto fc Maciell [l99p; ISeth et ail 120051; 



J0rgensenl |2000|; IWyse fc Gilmorel Il995l~ lHarris fc Harris! 



20021; iKoch et al.l 20061; ISaraiedini fc Jablonkal 1200a 
iMouhcine et all 120051 ; IWorthev et all 120051 ). Additionally, 
gas inflow over several Gyr is required to explain many 
characteristics of the SV, including stellar chemical com- 
positi ons and the pre s ent-d a y gas mass fraction a nd SFR 
(e.g., IChiappini et all l200ll ; IPortinari et~ai1 1 19981) . Other 
evidence for inflows includ es the kinematics of extra-p lanar 
HI gas in some spirals (|Fraternali fc Binnevl [2008), the 
low metallicities and large distances o f some MW HVCs 
l|Wakker et al.ll 19991 ; IRichter et al.ll200lT ). high velocity OVI 
absorption along vario us sight-lines through the MWs halo 
rtSembach et al. 20031) . and the prevalence of warps in HI 
discs l|Bosmalll99ir ). We refer the reader to ISancisi et alj 
(2008) for a review of observational evidence for gas inflow. 

There are theoretical expectations for gas inflows, as 
well. Cosmological simulations of galaxy formation in the 
ACDM framework predict disc galaxies to form from the ac- 
cretion of gas after the last major merger (lAbadi et al ]l2003l ; 
ISommer-Larsen et all 120031 ; iGovernato et al l |2004 12007V ). 
For galaxies with virialized dark halo masses < 10 12 M Q , the 
accreted gas is expected to be mostly cold (T ~ 10 4 ~ 5 K). 
On the other hand, upon entering the haloes of more mas- 
sive galaxies, most of the accreted gas experiences shock- 
heating to the virial temperature (T > lO 6 ^), creat- 
ing surrounding reservo irs of hot gas l|Keres et al. I 120051 ; 
iDekel fc Birnboiml [2006 ) . Thermal instabilities lead to the 
condensati on of cooler clouds, wh i ch then rain down on 
the discs (iMaller fc Bullock! |2004 ISommer-Larsen! 120061 ; 
iKaufmann et aLl2006l ). In some of the most recent N-body 
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simulations, these clouds have properties similar to the 
HVCs mentioned above l|Peek et al.l [2008). Furthermore, 
this idea is supported by the recent discovery of a hot 
(T ~ lO 6 /^) gas eous halo around the quiescent massive spi- 
ral, NGC 5746 (|Pedersen et alj|2006h . This gas is too hot 
to be heated by SNe in the disc and the disc SNe rate is 
too low to have created the reservoir through the outflow 
of disc gas. Finally, theoretical simulations suggest the need 
for extended accretion of dilute gas to keep discs from be- 
ing destroyed after a succession o f minor mergers with m ass 
ratios of 4:1 and even up to 10:1 |Bournaud et al.ll2007l ). 

Despite all this evidence, the precise nature and impor- 
tance of gas flows in the evolution of galaxies, particularly 
spirals, is still uncertain. To help improve the situation, in 
the present paper, we develop a chemophotometric CMD fit- 
ting method as an extension to the canonical method used in 
many other studies, but with the goal of e xamining the role 
of gas flows in M33's evolution. Following llkuta fc Arimotol 
(2002), we solve the chemical evolution equations to obtain a 
self-consistent SFH and AMR. We improve upon their work 
by allowing for gas inflow and outflow and by efficiently 
searching the full volume of parameter space to make a de- 
tailed and quantitative fit to the observed CMD. 

The photometric data we use in t his paper and the re- 
duction procedure were presented in iBarker et al.l |2007a, 
hereafter iPaper III) . In summary, three co-linear fields lo- 
cated in projection ~ 20 — 30' southeast of M33's nucleus 
were observed with the Advanced C amera for Surveys on 
board th e Hubble S pace Telescope. In IBarker et al.l (|2007bl . 
hereafter IPaper lit , we computed the SFHs for these fields 
using the canonical synthetic CMD fitting method with age 
and metallicity as free parameters. Because the outer two 
fields may have a no n-negligib le contribution from M33's 
halo or thick-disc fsee lPaper IlJ for a discussion), we restrict 
ourselves to analyzing the innermost field in the current 
study. This field has a projected galactic area of ~ 0.7 kpc 2 
and it li es at a deprojected radius of Rd p ~ 6 disk scale 
lengths (jFerguson et all 120071) or ~ 9 kpc, assuming a dis- 
tance of 867 kpc (iTiede et afll2004h. inclination of 56°, and 
position angle of 23° (|Corbelli fc Schneiderll997h . At this ra- 
dius in M33, the azimuthally averaged HI and stellar surface 
densities are, respectively, ~ 3 M^ pc~ 2 an d ~ 0.3 Mq pc -2 
dCorbefli fc Schneidenfl997l : ICorbellil[2003l ). 

In the next section, we outline the chemical evolution 
equations, which form the backbone of our models. We de- 
scribe in fj3]how we link these equations with the synthetic 
CMD fitting method to build a self-consistent model CMD. 
In SJH we present the results of fitting closed box and in- 
flow/outflow models to the data. We explore the effects of 
varying the model parameters in Sj5] In Sj6] we test the model 
predictions with independent observations and, in iJT] we 
compare the [a/Fe] vs. [Fe/H] relation in M33 with other 
LG systems. Finally, we summarize our results in SJ5] 



2 CHEMICAL EVOLUTION EQUATIONS 

To model the SFH of M33 in a self-consistent way we use 
the chemical evolution equations, which follow the varia- 
tion of the gas and stellar masses and the abundances of 
various elements. We use the instantaneous recycling ap- 
proximation (IRA) in which stellar lifetimes are treated as 



negligible compared to the age of the Universe for masses 
greater than the present-day main sequence turnoff mass 
(M > 0.8 Mq). This is a good approximation for follow- 
ing the gas and stellar masses and the abundances of el- 
ements like oxygen, whose primary producers are massive 
stars. However, the IRA breaks down for elements with sig- 
nificant contributions from SNe la because they are released 
to the ISM on longer time-scales. To follow the evolution 
of SNe la products we adopt the delayed produ ction ap- 
proximation (DPA) of lPagel fe Tautvaisienel (|l995l . hereafter 
PT95). In the DPA, there is a delay time, Td, between the 
birth of a stellar generation and the resulting SNe la explo- 
sions. 

Our models track the elements O, Mg, Si, Ca, Ti, and 
Fe. To minimize the effect of uncertainties in the yield 
of any one particular Q-element, we focus on their sum, 
which we refer to as a. Unless otherwise noted, [a/Fe] = 
[(0+Mg+Si+Ca+Ti)/(5 Fe)]. We do not follow carbon and 
nitrogen because they have a significant contribution from 
long lived low and intermediate mass stars, for which the 
IRA and DPA break down. In the future we plan to drop 
the IRA and DPA so that we can follow these two elements. 

We define the total mass surface density as E(t) = 
E s (£) + E g (t) where E s (t) and E 9 (t) are the stellar and 
gaseous surface densities, respectively, and t is the elapsed 
time of the model. We also define Ei (t) as the gas mass sur- 
face d ensity in t he for m of elemen t i. Combining the formal- 
ism of lTinslevI dl980h and |PT95l . the equations of chemical 
evolution under the IRA and DPA can be expressed as 

^ = fi(t)-fo(t) (1) 
2k = (l-RMt) (2) 
£± = -(l-R)1,(t) + f I (t)-fo{t) (3) 

= -(l-R)X i (t)i;(t)+y i (l-R)ij(t) 

+ y d ,i(i-R)i>(t-T d ) + fi,i{t)-foAt) ( 4 ) 

where tj) is the SFR, fi is the inflow rate, fo is the out- 
flow rate, and Xi = Ei/E 9 is the fraction of gas in the form 
of element i. The net stellar yield, y i; is the mass of newly 
synthesized element i instantaneously returned to the ISM 
by a stellar generation per unit mass locked up into stel- 
lar remnants. The delayed yield, yd,i, represents the delayed 
chemical yield of SNe la. The return fraction, R, is the frac- 
tion of mass returned to the ISM by a stellar generation. We 
calcul ate R from the yields of lvan den Hoek fc Groenewegenl 
d 19971) (for the case of a variable mass loss efficiency) and 
iPortinari et al.l dl998l ) . W e find that, ave r aged o ver all metal- 
licities, R = 0.235 for a iKroupa et ail d 19931 ) initial mass 
function (IMF). We note in passing that the precise value 
of R is inconsequential because the factor (1 — R) can be 
absorbed into the star formation efficiency parameter, e, de- 
scribed below. However, we explicitly include R to be more 
consistent with previous studies. 

Equs. [TJ — [4] represent conservation of total, gaseous, 
stellar, and elemental mass. The first term in Equ. [4] rep- 
resents the mass of element i that is originally present in 
the gas and lost to star formation, plus what is instanta- 
neously returned by the winds and explosive deaths of mas- 
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Table 1. Adopted chemical yields and solar abundances 



Element 


Vi 




Vi,d 








O 


6.68 x 10" 


3 


0.00 




9.546 x 10" 


3 


Mg 


4.53 x 10- 


4 


0.00 




5.144 x 10" 


4 


Si 


4.56 x 10" 


4 


7.82 x 10 


-5 


6.518 x 10" 


4 


Ca 


3.36 x HT 


-5 


1.08 x 10 


-5 


6.007 x 10- 


-5 


Ti 


1.49 x 10" 


-6 


2.56 x 10 


-7 


2.130 x 10- 


-6 


Fe 


3.23 x 10" 


4 


4.85 x 10 


-4 


1.155 x 10- 


-3 



a lAnders fe Grevessel Jl989h 

sive stars. The second term in Equ. [3] represents the instan- 
taneous restitution rate of newly synthesized element i while 
the third term is the delayed restitution of newly synthesized 
element i (from SNe la), which is zero for t < tj,. We adopt 
the semi-empirical yield s (show n in Table [TJ and time delay 
(r d = 1.3 Gyr) used by |PT95l . which they tuned to repro- 
duce the observed abundances and gas fractio n in the SV 
and then applied to the LMC and SMC (|PT9Sh . The fourth 
and fifth terms of Equ. [4] represent, respectively, the inflow 
and outflow rate of element i. The chemical composition of 
the inflowing gas is assumed primordial. The initial heavy 
element abundance and stellar mass are zero and, except for 
closed box models, the initial gas mass is zero, also. 

Motivated by the successes of c hemical evolution models 
appli ed to the MW's satellites (e.g.. lLanfranchi fe Matteuccil 
120041 '). we set the outflow rate proportional to the SFR by 
a constant factor, w, called the outflow efficiency. The com- 
position of the outflowing gas is the same as the ISM, which 
is assumed to be homogeneous and instantaneously mixed. 
In chemical evolution model s of LG dwarf galaxies the out- 
flow efficiency is ~ 1 — 10 (|Lanfranchi fe Matteuccil 120041 ; 
ICarigi et al.ll2006l ). so we restrict w to be less than 10. The 
idea behind this approach is that SNe inject kinetic energy 
into the ISM, possibly causing some gas to leave the sys- 
tem entirely. In reality, the outflow efficiency could depend 
on factors like the depth of the gravitational potential, the 
ejecta velocity, and the ISM's physical properties. A more 
detailed calculation including such properties is beyond the 
scope of the present study. We make no explicit distinction 
between radial or extra-planar gas flows. 

As is commonly done in chemical evolution models, 
we couple the SFR to the gas mass through the so-called 
Kennicutt- Schmidt (KS) relation, 



V»(t) = c52(t), 



(•») 



where e is the sta r formation effic iency. In the spirit of the 
pioneering work of lSchmidtl <|l959h . who found that star for- 
mation rate traced gas density in the MW. iKennicuttl (Il99&f ) 
measured mean gas masses and SFRs within the optical 
radii of ~ 60 normal spirals and within the central regions 
of ~ 30 starburst galaxies and found e = 0.25 ± 0.07 and 
re — 1.40 ± 0.15 0- Looking at the two subsamples individ- 
ually, he found k = 2.47 ± 0.39 for the normal spirals and 
k — 1.40 ± 0.13 for the starbursts. 

This global star formation relation can be compared to 



1 on a scale where Mq pc 2 Gyr 1 are the units of ip and 
Mq pc -2 are the units of S 9 



a local one which refers to gas densities and SFRs at individ- 
ual points or in azimuthally averaged bins within a galaxy. 
Several studies of nearby sys tems have found that, when ex- 
amined locally , re sa 1 .2 — 3.5 (Gottesman fe Weliachewlll977l; 
Wong fc Bhtzl |2002|; iKennicutt et. al.1 12007| ; iBoissier et al.1 
20031 ; iMisiriotis et al.ri2006h . The cause of the observed 



variation in k is unclear but possibilities include non- 
axisymmetric profiles or uncertainties in the extinction cor- 
rection, the CO-H2 conversion factor, flux calibration, and 
the conversion to SFR. The variations may be due to intrin- 
sic properties of the galaxies, as there is some evidence that 
molecule-rich galaxies, typically massive spirals and star- 
bursts, have lower re than molecule-poor galaxies, typically 
low surface brightness and dwarf systems. Interestingly, M33 
falls into this latter cate gory with a total molecular fraction 
of ~ 0.1 |Corberlill2003h . 

The dependence of re on molecular fraction is not too 



ular j 


?as quite well (e.g., 


Murgia et al. 


2002 


; Hever et al. 


2004; 


Matthews et al]|2005 


; iLerov et al.l 


2005; 


Gardan et al. 


200?t). Indeed, the studies mentioned above also found that. 



when considering the molecular gas alone, re ~ 1.4, with 
much less galaxy-to-galaxy variation than when considering 
the total or atomic gas. It is beyond the scope of the present 
study to track the time evolution of the molecular gas frac- 
tion, but this is one possible avenue for future improvement 
to the techniques presented here. 

With these considerations in mind, an appropriate em- 
pirical KS relation to use wo uld be one de r ived f or M33 in- 
volving the total gas density. iHever et all |2004) measured 
k = 3.3 ± 0.07 and e = 0.0035 ± 0.066 using the infrared 
luminosit y to estimate M33's S FR profile inside Rd p ~ 30'. 
Similarly, IBoissier et. all (|2007h studied M33's KS relation 
by measuring its UV surface brightness and t ransla ting that 
to a SFR profile. Although IBoissier et. al.l (|2007l ) did not 
report specific values, w e estimate by eye fr om their figures 
a similar value for re as IHever et all (|2004n , but a smaller 
value for e by a factor of ~ 3. Adopting the IHever et al.l 
(2004) values for e and re leads to very poor fits of the CMD 
because the distributions of stellar age and metallicity are 
skewed too high and low, respectively (see fj4]2}. Further 
tests show that our models are more sensitive to reasonable 
changes in e than in re, so we adopt re = 3.3 for all models 
and allow e to be a free parameter (constant in time) re- 
stricted to the range 0.0001 — 10.0. We note that adopting 
re = 1.4 does not significantly change our general conclu- 
sions because the gas mass surface density ~ 1 Mq pc -2 
throughout much of the system's evolution. Finally, we do 
not inclu de a star formation thresh old (see the discussion in 
IPaper ij and IBoissier et. all (|2007l ) for reasons why). 

Our chemical models can be summarized as follows. Gas 
inflow deposits gas into the system and drives star formation 
via the KS relation. Because the inflowing gas has primor- 
dial composition, the ISM metallicity is always below what 
it would be without inflow. About 25% of the gas that goes 
into making each stellar generation is instantaneously re- 
turned to the ISM by stellar winds and SNe II, while ~ 75% 
is locked up into stellar remnants like white dwarfs, neu- 
tron stars, and black holes. The star formation enriches the 
ISM with metals, but also drives an outflow of gas from 
the system. All else being equal, the presence of a gas out- 
flow quenches the SFR, slows the chemical enrichment, and 
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Table 2. Fit qualities, distances, and extinctions of the Padova solutions 



Model Q xl v {m — M) a Ay <j 



Closed box 


29.13 


2.42 


1743 


24.50 


0.05 


0.25 


0.03 


Ex. inflow 


19.28 


2.03 


1741 


24.56 


0.07 


0.21 


0.05 


Sand, inflow 


15.37 


1.99 


1741 


24.60 


0.05 


0.18 


0.04 


Double Ex. inflow 


9.16 


1.81 


1739 


24.60 


0.05 


0.18 


0.04 


Truncated inflow 


8.68 


1.78 


1739 


24.60 


0.05 


0.20 


0.06 


Free inflow 


8.49 


1.77 


1739 


24.60 


0.05 


0.20 


0.06 



Table 3. Star formation and outflow efficiencies of the Padova solutions 



Model 


log(e) 


Chi 


°"lo 


w 


Chi 


Clo 


Closed box 


-2.29 


0.11 


0.05 








Ex. inflow 


-1.43 


0.72 


0.44 


0.66 


1.11 


0.52 


Sand, inflow 


-0.69 


0.76 


0.64 


0.99 


0.67 


0.52 


Double Ex. inflow 


0.77 


0.48 


1.28 


1.44 


0.47 


0.80 


Truncated inflow 


-0.47 


0.55 


0.58 


0.98 


0.55 


0.56 


Free inflow 


-0.21 


1.05 


0.65 


1.21 


0.75 


0.78 



suppresses the metallicity below what it would be without 
outflow. The ISM metallicity generally increases with time 
due to stellar nucleosynthesis, but a short, rapid increase 
in the inflow rate (IFR) can decrease the ISM metallicity 
temporarily. In the absence of all gas flows, an initial gas 
reservoir must be present to start SF, which continuously 
declines thereafter. 



3 METHOD 

We start with a set of model parameters (e.g., e, k, w, 
and //), which we use to solve Equs. [T] — [3] We integrate 
these equations using a 4th-order Runge-Kutta method with 
a time-step of 2 x 10 6 yr. The numerical integration was 
checked against several analytic solutions and the typical 
fractional error was ~ 10~ 6 . The resulting chemical evolu- 
tion model specifies the total and gas mass surface densities, 
SFR, and abundances of the various elements as functions 
of time. 

The age-metallicity plane is divided into logarithmic 
bins of width 0.25 dex in age and 0.3 dex in metal abun- 
dance. To c over this p lane, we use the same set of synthetic 
CMDs as in lPaper Hit each of which represents the predicted 
photometric distribution of stars in the corresponding age 
and metallicity bins. Ea ch of these CMDs was cre ated with 
the IAC-star program (|Aparicio fe GallartJ liooi) . We cal- 
culate the total stellar mass formed in each age bin from 
the chemical model and the corresponding mass-weighted 
mean global metallicity from Equ. [6] described below. If the 
metallicity falls outside the CMD library, it is capped at 
the appropriate limit. In each age bin, the total stellar mass 
is split up into the two synthetic CMDs which bracket the 
mean metallicity. More weight is given to the CMD whose 
central metallicity is closer to the mean metallicity. The fi- 
nal model CMD is the superposition of all synthetic CMDs 
weighted by their amplitudes. The data and model CMDs 
are divided into square bins 0.1 mag on a side and compared 
using the maximum likelihood ratio, T, for a Poisson distri- 



bution. The model parameters are changed and the whole 
process repeats until T is minimized. 

The quality of the fit is formally measured by the pa- 
rameter, Q, which gives the difference between T and its 
expectation value in units of its standard deviation. A Q 
value of 0.0 indicates a perfect fit, while, for exampl e, a value 
of 2.0 indicates a 2a departure from a perfect fit (|Dolpmnl 
l2002l . |Paper"lll . Because many readers may be more famil- 
iar with x 2 tha n Q, we al so calculate the reduced \ 2 of the 
best fit. As in iPaper IlJ, we simultaneously solve for dis- 
tance and extinction with a fixed grid spanning the ranges 
(m - M) = 24.50 - 24.80 and A v = 0.10 - 0.25 in steps of 
0.10 and 0.05 mag, respectively. The global best fit (referred 
to simply as the best fit) is the weighted average of those 
distance/extinction combinations whose solutions lie within 
la of the best individual solution. 

W e use the program, StarFISH l|Harris fc Zaritskvl 
l200lT). after incorpo rating the genetic algorithm, PIKAIA 
(|Charbonneaul 11995). to efficiently search the full volume 
of parameter space. Briefly, this algorithm randomly gen- 
erates an initial population of solutions, which is evolved 
through successive generations under the action of natu- 
ral selection, breeding, inheritance, and random mutation 
to find the global solution. We ran PIKAIA for 200 genera- 
tions in the steady-state-replace-random reproduction mode 
with creep mutation enabled and 50 individuals in each gen- 
eration. Then, the downhill simplex routine of StarFISH 
was started from the best-fitting solution found by PIKAIA. 
This hybrid approach proved more efficient at locating the 
global solution than either method individually. 

Because we are no longer solving for the amplitudes of 
the basis populations, we had to modify the way StarFISH 
calculates the confidence intervals. For each parameter, we 
take small steps in the positive and negative directions away 
from its optimum value. After each step, we allow the down- 
hill simplex to re-converge while holding that parameter 
fixed and allowing the others to vary. The process is re- 
peated until the lcr limit of the fitting statistic was reached. 

To estimate the systematic errors in the stellar evolu- 
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Table 4. Fit qualities, distances, and extinctions of the Teramo solutions 



Model Q xl v {m — M) a Ay <j 



Closed box 


33.26 


2.41 


1743 


24.63 


0.08 


0.17 


0.04 


Ex. inflow 


22.61 


2.18 


1741 


24.60 


0.05 


0.10 


0.03 


Sand, inflow 


15.08 


2.02 


1741 


24.70 


0.05 


0.10 


0.03 


Double Ex. inflow 


9.26 


1.90 


1739 


24.70 


0.05 


0.10 


0.03 


Truncated inflow 


9.59 


1.90 


1739 


24.70 


0.05 


0.10 


0.03 


Free inflow 


8.93 


1.89 


1739 


24.70 


0.05 


0.10 


0.03 



Table 5. Star formation and outflow efficiencies of the Teramo solutions 



Model 


log(e) 


Chi 


°"lo 


w 


Chi 


Clo 


Closed box 


-1.75 


0.16 


0.15 








Ex. inflow 


-0.19 


0.70 


0.28 


0.15 


0.72 


0.15 


Sand, inflow 


0.33 


0.64 


0.84 


0.82 


0.18 


0.57 


Double Ex. inflow 


0.99 


0.00 


0.77 


0.78 


0.14 


0.30 


Truncated inflow 


0.30 


0.70 


0.31 


0.40 


0.52 


0.25 


Free inflow 


1.00 


0.00 


0.57 


0.77 


0.16 


0.20 



tionary tracks we cre a ted tw o set s of synthetic CMDs us - 
ing the iGirardi et ail (|2000h and iPietrinferni et all (12004m 
tracks, which we respectively designate as Padova and Ter- 
amo (also referred to as BaSTI in the literature). The con- 
version from the theor etical to the observation al plane is 
accomplished with the ICastelli fc Kuru cz (2003) library of 
bolometric corrections. 

At this point it is important to note that these stel- 
lar tracks have scaled-solar abundances, but the elemental 
abundance distribution in forming stars changes with time. 
This is particularly important for the RGB, whose temper- 
ature is determined mostly by the abundances of elements 
with low first ionization potentials (Fe, Mg, Si, S, and Ca), 
which are the primary atmospheric electron donors. Ide- 
ally, the stellar tracks would exist for all possible elemental 
abundance distributions and we would pick the appropriate 
tracks at each time step of the chemical model. However, 
the computational effort required to generate such a com- 
prehensive set of tracks is too prohibitive. Only recently have 
groups begun to compute tracks for many different masses 
and metallicities and different levels of a-enhancement, as- 
signing similar enhancement to al l the a-elements, but on 
a very coarse grid in [a/Fe] (e.g., VandenBerg et all 120061; 
IPietrinferni et al.ll2006l ; IDotter et al.ll2007r ). 

ISalaris et all (| 19931 ) found that a-enhanced tracks and 
isochrones are well reproduced by scaled-solar ones with 
the same global metallicity provided the enhancements 
are similar for all the a-elements. Subsequent studies 
demonstrated that this result begins to break down for 
Z ;> 0.003 (e.g. ISalaris fc Weissl ll99S| ; ISalasnich etail 



the fo rmalism of ISalaris et alj (| 1993T ) and IPiersanti et al.l 
< |2007h . 



(6) 



l2000l I VandenBerg et al.l l2000l : iKim et al.l l2002h . We have 
checked its applicabi lity using several of th e most recent 
isoch r one databases |Pietrinferni et al. 2006; P otter et al.l 
l2007l ; IVandenBerg et all I2006T I and we find that over the 
range of ages, [Fe/H], and [a/Fe] most appropriate for our 
data, a scaled-solar isochrone in the I vs. (V — I) plane is 
within ~ 0.1 mag of an a-enhanced isochrone with the same 
global metallicity. Therefore, we account for a-element en- 
hancements by calculating the global metallicity following 



[M/H] = [Fe/H] +log(a a 10 [Q/Fcl +b Q ) 

In Equ. HO a a = J2i( X i/ Z )& ~ °- 7 for 1 = °> Ne > M S> 
Si, S, Ca, and Ti and b a = 1 — a Q . Note that this relation 
implicitly assumes no enhancements in elements other than 
the a-elements. This is not a significant problem, however, 
since the a-elements comprise most of the metals by mass. 
In any case, our conclusions are insensitive to the precise 
value of the correction factor to [Fe/H] in Equ. HJ] because it 
generally amounts only to < 0.2 dex. 

In principle, the bolometric corrections, color transfor- 
mations, and evolutionary lifetimes also depend on [a/Fe], 
but in practice, our data are not signif icantly affected by 
these dependencies. ICassisi et alj (|2004T ) showed that the 
bolometric corrections and color transforma tions depend 
neglig ibly on [a/Fe] in V and redder bands. IDotter et al.l 
(2007) investigated the effect of abundance variations on 
stellar evolutionary models and found that when [a/Fe] = 
0.3, the MS lifetime is decreased by ~ 5%. We have com- 
pared the scaled-sola r and a-enhanced Teram o tracks, which 
have [a/Fe] ~ 0.4 (|Pietrinferni et all I2006T ). and we find 
that, in general, the evolutionary phase lifetimes differ by 
< 10%. These effects are likely to be even smaller for our 
data since we derive —0.2 < [a/Fe] < 0.2 for almost all ages. 



4 RESULTS 

We tested several different inflow/outflow scenarios using 
both sets of stellar tracks. In general, they give similar re- 
sults, so we only show those of the Padova tracks. In the 
text, we describe any significant differences, since they can 
help us gauge the systematic errors due to the tracks them- 
selves. The top row in Figs. [T]- [2] shows the residual CMD 
of a particular scenario on a scale where white and black 
correspond, respectively, to an excess and deficit of model 
stars at the 3<r level. The bottom row in these figures shows 
the best-fitting SFH as the upper line and inflow history 
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Figure 1. Star formation history and inflow history results using the Padova tracks. The top panel in each column shows the residual 
CMD on a range of ±3cr where positive (white) residuals indicate an excess of model stars. The bottom panel in each column shows the 
star formation history as the upper line and inflow history (scaled by a factor of 0.1 and integrated over the total field area) as the lower 
line. From left to right, the col umns corr espond to the closed box, exponential inflow, and Sandage inflow models, respectively. The last 
column shows the results from iPaper Ilfl, for which age and metallicity were free parameters. 



(IFH) divided by 10 as the lower line. The IFH is integrated 
over the total field area. Th e last column in these figures 
gives the result of IPaper IlJ, in which age and metallicity 
were free parameters. Tables [5] and [3] give the the fit quality 
(Q), reduced \ 2 (x?)i number of degrees of freedom (v), and 
mean distance modulus and extinction with their respective 
la uncertainties. Tables [3] and [S] give the mean star forma- 
tion and outflow efficiencies and their upper and lower la 
uncertainties. The mean values reported in the tables are av- 
erages of the acceptable solutions in the distance/extinction 
grid as explained in §3. 

4.1 Closed Box Models 

We began by testing the canonical closed box model in which 
the inflow and outflow rates were identically zero. The sys- 
tem was initially composed entirely of gas. The total mass 
remained constant throughout its evolution while the gas 
mass and SFR decreased monotonically with time. The best- 
fitting closed box model is displayed in the first column of 

Fig.m 

There are several significant discrepancies between the 
model and data CMDs. First, the model predicts too much 
star formation at ages < 1 Gyr which causes an overabun- 
dance of model main sequence (MS) stars on the blue plume 
at (V — I) ~ 0.0. Second, the model has too few stars in a 
region below the blue plume centred at (V — I) ~ 0.5 and 



/ ~ 26.5. In the global solutions of IPaper Hit this region is 
dominated by stars with ages ~ 2 — 8 Gyr, indicating the 
closed box model has too little star formation at these ages. 
Also, the RGB of the model is too blue, indicating a mean 
metallicity that is too low. 

The Teramo model exhibits similar discrepancies, but 
the RGB appears slightly too wide, possibly indicating an 
excessively large metallicity spread, and it has too many 
stars overall. The Teramo model also has too many stars on 
the blue horizontal branch, which signals that there is too 
much star formation at the oldest ages. It cannot be due to 
the metallicity of the oldest bin being too low because it is 
already close to that of the IPaper IlH solution, which does 
not show such a discrepancy. 

The fit qualities of the closed box solutions are > 20<r 
worse than the solutions of IPaper lit which had Q values 
of 6.64 (Padova) and 6.02 (Teramo). This is not surprising 
since our closed box model has only two free parameters, 
namely, the star formation efficiency and the total mass. 
However, as we will see below, the addition of gas inflow 
and outflow significantly improves the fits with only 2 — 4 
more free parameters, indicating this region in M33 proba- 
bly did not evolve as a closed box. Because of the reduced 
number of free parameters, the confidence intervals of the 
clos ed box so lutions are significantly smaller than those of 
the IPaper fill solutions. This demonstrates that simultane- 
ous CMD and chemical evolution modelling can be used to 
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alleviate the age-metallicity degeneracy inherent in broad- 
band stellar colors. 



4.2 Inflow and Outflow Models 

As described in iJU there is evidence that galaxies do not 
evolve as closed boxes. An exponential inflow rate is one of 
the simplest and most common forms used in the literature, 
so it is instructive to see how well it can explain M33's SFH. 
Accordingly, we solved for the inflow time-scale and the ini- 
tial inflow rate. We also included gas outflow by allowing a 
nonzero outflow efficiency, w. The second column in Fig. [T] 
shows the exponential inflow model for the Padova tracks. 

This new model provides a better fit to the data than 
the closed box model, but it still exhibits some large dis- 
crepancies with the data. In fact, these discrepancies are 
very similar to those of the closed box model, but the mag- 
nitude of the residuals has been lessened. There is still too 
much star formation at ages < 1 Gyr and too little in the 
range 2 — 8 Gyr. The fit quality is ~ 15a worse than the 
IPaper IIP solutions. 

Some numerical simulations of structure formation 
within the ACDM framework predict that the average 
mass accretion rates of dark matter haloes are initially 
small, grow to some maximum, and decline thereafter (e.g., 
van den Boschl |2002| ; IWechsler et ail I2002T ). With this in 
mind, we also investigated another function for the IFR 
which has a delayed maximum, fi(t) oc t exp(—t 2 /2r 2 ), 
where r is the time between when the inflow starts and when 
it peaks. This fu nction, which w e refer to as Sandage inflow, 
was first used bv lSandaed (|f 9861 ) to describe the variation in 
SF H with galaxy mo r pholo gy and later explicitly presented 
by iMacArthur et al.l l |2004 ). By shifting the bulk of the in- 
flow toward younger ages (i.e., earlier lookback times), it 
allows for more star formation at correspondingly younger 
ages. The result is plotted in the third column of Fig. Q] 
The morphologies of the residuals are similar to the expo- 
nential inflow model, but their m agnitudes a re smaller. The 
fit quality is ~ 9a worse than the lPaper flj solutions. 

The Teramo models show a qualitatively similar behav- 
ior - the exponential and Sandage functions do a progres- 
sively better job at reproducing the observed CMD, but still 
do not get the overall age distribution correct. The main 
drawback of these functions is that the IFR today cannot 
easily be varied independently of the IFR at intermediate 
ages (~ 2 — 8 Gyr). These functions just do not provide 
enough freedom to describe the true IFH which may not 
be adequately characterized by just two parameters or a 
smooth function. 

These results led us to try three less restrictive inflow 
models. The first was a double exponential model described 
by four parameters: a growing time-scale, a decaying time- 
scale, a transition time between the growing and decaying 
modes, and the IFR at the transition time. The second was 
a truncated model described by four parameters: the initial 
IFR, the IFR at a model time of 7 Gyr, a truncation time 
when the inflow ends, and the IFR at the truncation time. 
In the third model, which we called free inflow, we approxi- 
mated the true IFH with a discrete function by dividing the 
entire age range into 4 bins each with a possibly different 
but constant IFR. We experimented with different binning 
schemes and settled upon a compromise between the de- 



sire to have a reasonable computing time and small error 
bars (requiring fewer bins) and the competing desire to pro- 
duce complicated SFHs (requiring more bins). The free in- 
flow model allows discontinuous jumps in the IFR, but these 
should not be interpreted too literally. This is the same ap- 
proach that we take in approximating the true SFH with 
9 logarithmic bins, each with constant SFR. Internal tests 
have shown that these inflow functions can reasonably cap- 
ture broad trends in the true IFH and chemical evolution, 
especially at ages < 7 Gyr. 

The solutions using these three inflow models are dis- 
played in Fig. [5] The free inflow model provides the highest 
quality fits. However, the double exponential and truncated 
inflow models are < la worse and exhibit qualitatively sim- 
ilar results. The discrepancies between model and data ex- 
hibited previously with the exponential and Sandage inflow 
models are almost completely erased. Most of the discrep- 
ancies that do remain are similar to those exhibited by the 
solutions in IPaper Hit leading us to conclude that they are 
mostly caused by inaccuracies in the stellar tracks. However, 
the fit qualities are still ~ 2 — 4<r worse than the IPaper IlH 
solutions. This difference could arise from the approxima- 
tions made in producing the chemical evolution models and 
the uncertainties inherent in the stellar yields. Second, age 
and metallicity are no longer completely free parameters so 
errors in the stellar tracks cannot be as easily hidden by 
arbitrary combinations of age and metallicity. Third, we 
do not account for a systematic change in the metallicity 
spread with age. Finally, we have not allowed for an inflow 
or outflow of stars with other ages and metallicities into 
the single coherently evolving region we are modelling. This 
mixing of stars, perhaps caused by gravitational interactions 
with spiral arms, could artificially increase the spread in the 
age-metallicity relation and it could mean the stars in this 
region originated, on ave r age, at another lo cation in M33 
l|Sellwood fc Binnevll2002] : iRoskar et alll2008ft . 

In Fig. [3] we summarize the distribution of stellar ages 
and metallicities in the three best inflow models for the 
Padova tracks. The panels show the (a) SFH, (b) age cumu- 
lative distribution function (age CDF), (c) Z AMR, and (d) 
Z metallicity distribution function (Z MDF) of all stars ever 
formed. The dotted and dashed lines are, respectively, the 
double exponential and truncated inflow models. The solid 
line with diamonds and error bars corresponds to the free in- 
flow model. The open circles are the results from IPaper Ilj . 
Note that only the free inflow model errors are shown for 
clarity. Additionally, we show the IFH and inflow CDF of 
these three models in Fig. 2] 

In the best three inflow models, ~ 50 — 60% of the gas 
accretion takes place between 3 and 7 Gyr ago. Anything 
less would produce too few intermediate-age stars. Second, 
at most ~ 10% of the gas was accreted in the last 3 Gyr. 
Any amount in excess of that would lead to a recent SFR 
that is too high and produce too many young stars on the 
blue plume MS at (V - I) ~ 0.0. Third, the preferred out- 
flow efficiency is ~ I ± 0.5, which is smaller than the typ- 
ical values of ~ 5 — 10 estimated for dwarf galaxies i n the 
LG l|Lanfranchi fc Matteuccil liooi : ICarigi et alj|2006h . Fi- 
nally, the results f or the SFH , age CDF, AMR, and Z CDF 
are similar to the IPaper IlH so lutions. This lends support 
to our conclusions in IPaper Hit and suggests they are not 
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Figure 2. Same as Fig.[T] but the first three columns show the double exponential, truncated, and free inflow models, respectively. 



significantly affected by unphysical combinations of age and 
metallicity. 

The preferred values for e tend to be of order ~ 0.1 — 1.0. 
The resulting gas depletion timescale, r 9 = T, g /T,sFR ~ 
3 — 10 Gyr, is on the low end of the equivalent timescale in 
dwar f galaxies, which ranges from a few to several tens of 
Gvr dTavlor fc Webster1l2005l ; iKarachentsev fe Kaisinll2007l; 



Calura et al. 2008). An e value as low as 0.0035 l|Hever et al.l 
20041 ) is strongly disfavored because, as Fig. ^demonstrates, 
the resulting distributions of stellar age and metallicit y are 
skewed to higher and lower values than the IPaper IlH solu- 
tion. If, instead, we fix e to have a higher value, like 0.084, as 
implied by the local current gas density and SFR in our field, 
then the age and metallicity distributions are more reason- 
able. T his value of e is only 1.2a larger than the lHever et alj 
(2004) value and only 1.1a lower than the best-fitting value 
of 0.62 in the free inflow model. The variation in e among the 
best three inflow models indicates the exact parametrization 
of the IFH can affect it by almost an order of magnitude. 
Based on the tests in fj5j a similar uncertainty in e arises 
from various assumptions built into the chemical models, 
like the stellar yields and the metallicity of the inflowing 
gas. 

By incorporating the chemical evolution equations into 
the CMD fitting, we can extract more information from 
the solutions and make more predictions that can be tested 
against observations. Fig.[6]summarizes some predictions for 
the three best inflow models. Each figure shows (a) the gas 
mass averaged over each of the 9 SFH bins, (b) the stellar 
mass- weighted ([a/Fe]) and {[Fe/H]) of all stars formed in 
each SFH bin, (c) the Fe AMR, and (d) the Fe MDF. In 



panel (b), the horizontal and vertical lines denote {[a/Fe]) 
and ([Fe/H]) of all stars ever formed, respectively. The line 
types are the same as in Fig. [3] 

Because of the presence of gas inflow, the gas mass rises 
in the early evolutionary stages and reaches a maximum 
several Gyr later. The gas mass begins to decline as the 
inflow rate becomes small but star formation and outflow 
continue. Interestingly, the total mass actually decreases in a 
few cases when the outflow rate exceeds the inflow rate. The 
mean [a/Fe] of the oldest age bin is ~ 0.2 ±0.1 while that 
of the youngest bin is ~ —0.1 ±0.1. Because the majority of 
stellar mass formed within the oldest three age bins, ([a/Fe]) 
of all stars ever formed is ~ 0.1, even though the majority 
of age bins have ([a/Fe]) < 0.0. Note that the [a/Fe] vs. 
[Fe/H] relation is not always single- valued, since [a/Fe] can 
increase and [Fe/H] can decrease with time depending on 
the precise interplay between gas flows, the SFR, and the 
SN la explosion rate (~see l7.ll) . 

The Teramo models show an overall similar history as 
the Padova models, with ~ 50% of the total inflow tak- 
ing place in the last 7 Gyr and < 5% in the last 3 Gyr. 
However, the inflow hiatus present in the Padova free inflow 
model between ~ 6 and 9 Gyr is not present in the Teramo 
model. The mean [a/Fe] is also similar between the two sets 
of tracks. The Teramo models predict gas masses smaller by 
a factor of ~ 2 and metallicities ~ 0.2 dex higher at all ages. 
This meta llicity difference is somewhat larger here than in 
IPaper IlJ . We believe this is mostly due to the metallicity 
of the 6.2 Gyr age bin in th e Teramo s olution now being 
~ 0.3 dex larger than it is in IPaper Hit The overall faster 
enrichment of the Teramo solutions compared to the Padova 
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Figure 3. (a) Star formation history, (b) age cumulative distribution function, (c) age-metallicity relation, and (d) metallicity distribution 
function of the acceptable inflow solutions using the Padova tra cks. Each panel shows the double exponential (dotted), truncated (dashed), 
and free inflow (solid) models and the results from lPaper IlJ (open circles). Only the free inflow model errors are shown for the sake of 
clarity. 



solutions requires higher fitted values of e. In a couple of the 
inflow models, e reaches its upper limit of 10. We repeated 
these fits with an upper limit of 100 and the resulting fit 
qualities were improved by ~ 1 — 2a, but the details of the 
SFH, IFH, and chemical evolution were not significantly dif- 
ferent. 



5 VARYING THE MODEL PARAMETERS 

There were several parameters in our chemical models which 
we fixed to match observations or be consistent with previ- 
ous studies. These parameters included the initial gas mass, 
initial chemical composition, composition of the inflowing 



gas, the KS relation exponent, and the stellar yields. To es- 
timate their effects on our results, we varied each of these 
parameters and repeated the free inflow fit. This resulted in 
several new models which we refer to as A through K and 
which had the following parameter changes: (A) an initial 
unenriched gas reservoir with E s (0) = 1.0 M Q pc -2 , (B) the 
same as A but the gas reservoir was pre-enriched to [Fe/H] 
= —1.3, (C) the inflowing gas was pre-enriched by SNe II 
to [Fe/H] = —1.3, (D) a superposition of A and C, (E) a 
superposition of B and C, (F) yi was increased by a factor 
of 2, (G) yi was decreased by a factor of 2, (H) k — 1.4, 
and (I) a different binning scheme was used. We also ran 
two additional tests in which e was allowed to vary with 
time: (J) a closed box model with an initial metallicity of 
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Figure 4. (a) Inflow history and (b) inflow cumulative distribution function of double exponential (dotted line), truncated (dashed line), 
and free inflow models (solid line) using the Padova tracks. For clarity, only the free inflow model errors are shown. 
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Figure 5. (a) Star formation history and (b) age-metallicity relation using two independently estimated values for M33's star formation 
efficiency, e = 0.0035 (dotted line) and 0.084 (dashed line; see text for details). The original free inflow model (sol id line with diamonds) 
treated e as a free parameter, resulting in a best- fitting value of 0.62. Open circles are the results from lPaper IlJ . Only the original free 
inflow model errors are shown for the sake of clarity. 



[Fe/H] = —1.3 and (K) a free inflow model with a constant 
IFR (i.e., one inflow bin rather than 4). While these tests 
are by no means exhaustive, they give a rough sense of the 
potential systematic errors that could be introduced by the 
assumptions we made in our chemical models. 

Figs. [7] -[9] compare the original best- fitting free inflow 
model (line with diamonds and error bars) to the new models 
(lines without error bars). The open circles are the results 



from lPaper IlJ . Note that only the original free inflow model 
errors are shown for clarity. 

In general, the new results are not significantly different 
from the original results and the fit qualities are unchanged 
to within ~ la. The SFH, age CDF, MDF, and Z CDF are 
the least affected by the n ew parame ter values and they re- 
main close to the results of IPaper IlJ . More variation occurs 
in the IFH, inflow CDF, gas mass evolution, and [a/Fe] vs. 
[Fe/H] relation. The IFR and gas mass in any age bin can 
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Figure 6. Gas mass and chemical composition in the double exponential (dotted line), truncated (dashed line), and free inflow models 
(solid line) using the Padova tracks, (a) Evolution of the gas mass surface density, (b) [ct/Fe] vs. [Fe/H] relation, (c) age-metallicity 
relation for iron, and (d) metallicity distribution function for iron. Only the free inflow model errors are shown for the sake of clarity. 
The horizontal and vertical line segments in panel (b) denote, respectively, the mean [a/Fe] and [Fe/H] of all stars ever formed. 



show variations up to a factor of ~ 2. The presence of an 
initial gas reservoir, as in models A, B, D, and E, decreases 
the amount of inflow required in the oldest bin because oth- 
erwise the resulting SFR would be too high. This means 
that, all else being equal, the fraction of total inflow occur- 
ring at younger ages (< 7 Gyr) in the original models is 
a lower limit. Increasing [Fe/H] of the initial gas reservoir 
only increases ([Fe/H]) and ([a/Fe]) of the oldest age bin. All 
else being equal, increasing [Fe/H] of the inflowing gas shifts 
([a/Fe]) of all age bins upward ~ 0.1 dex and shifts the ini- 
tial and final ([Fe/H]) upward and downward by ~ 0.2 dex, 
respectively. Models C - F have some of the largest IFRs 
and gas mass densities. In model F, a larger IFR of pri- 
mordial gas is needed to balance the increased stellar yields. 



When the inflowing gas is not primordial, as in models C - 
E, more of it is required to counteract the ISM enrichment 
due to stellar nucleosynthesis. 

One of the largest sources of uncertainty in chemical 
evolution models is the stellar yields. Changing the instan- 
taneous yields, j/i, by a factor of 2, as in models F and G, 
shifts the [a/Fe] vs. [Fe/H] relation vertically by ~ 0.1 — 0.2 
dex with only a small change to the relation's tilt. Increasing 
or decreasing the delayed yields, yi,d, has the same effect on 
the [a/Fe] vs. [Fe/H] relation as decreasing or increasing, re- 
spectively, the instantaneous yields. This fact arises because 
the a-elements contribute the majority of yi while iron con- 
tributes the majority of y iy d- Decreasing yi or increasing j/; it j 
appear to be the only ways to significantly move the en- 
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Figure 7. Effects of varying the model parameters on the star formation history and age-metallicity relation using the Padova tracks. 
The solid line with diamo nds is the original free inflow model and the other lines are tests A through K (sec text for details). Open 
circles are the results from lPaper IlJ . Only the original free inflow model errors are shown for the sake of clarity. The panels are identical 
to Fig. GS 



tire [cv/Fe] vs. [Fe/H] relation down relative to the original 
model. 

To demonstrate the binning effects on the results, Model 
I has a different binning scheme from the original free inflow 
solution. Instead of 4 bins, Model I uses the same 9 loga- 
rithmically spaced bins as the SFH. This model is within 
la of the original free inflow model and has a similar inflow 
history and chemical evolution to within the errors. The two 
most important differences are that, first, the inflow hiatus 
appearing in the original model is no longer present. This 
hiatus was possible in the original model, and not the new 
one, because of the overlap between the 2nd oldest SFH bin 
and 3rd oldest IFH bin. Second, the new binning scheme al- 
lows a large inflow burst in the youngest two bins (~ 80 — 250 



Myr) amounting to ~ 10% of the total inflow. During this 
burst, the IFR increases by a factor > 1000 over the pre- 
vious two bins. This burst occurs in order to repr oduce the 
apparent SFR burst in the youngest bin seen in thelPaper IlJ 
solution. As the tests and discussion in iPaper IlJ bear out, 
the reality of such a SFH burst is highly suspect because 
the SFR amplitudes in the youngest couple of age bins are 
susceptible to low number statistics. 

In all the preceding results, we have attributed varia- 
tions in SFR primarily to variations in IFR because the star 
formation efficiency, e, was constant with time. A natural 
question to ask is, can the SFR variations be due, instead, 
to variations in e? Allowing e to vary weakens the link be- 
tween SFH and IFH, but does not completely decouple them 
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Figure 8. Effects of varying the model parameters on the inflow history using the Padova tracks. The solid line with diamonds is the 
original free inflow model and the other lines are tests A through K (see text for details). Only the original free inflow model errors are 
shown for the sake of clarity. The panels are identical to Fig. [4] 



because e and the IFR have opposite effects on the ISM 
metallicity. A sudden increase in the IFR dilutes the metals 
in the ISM, decreasing the enrichment rate, whereas a sud- 
den increase in e increases the SFR and, hence, enrichment 
rate. Phenomena like gravitational perturbations and stellar 
feedback could potentially change e, which may more gener- 
ally evolve according to local gas phase properties like tem- 
perature, density, pressure, chemical composition, molecu- 
lar fraction, and interstellar radiation field strength (e.g . 
Elmegreenlll993l; ISommer-Larsen et al. 20031; ISchavei 12004; 



Schave fc Dalla Vecchiall200Sl ; lRobertson fc Kravtsovf2008h 



Models J and K respectively demonstrate that a closed 
box or constant IFR model provide an acceptable fit to the 
CMD if we allow e to vary with time. However, the amount 
of variation required (an r.m.s. deviation over time ~ 0.5 
dex) may be larger than the intrinsic scatter in the KS rela- 
tion of other spirals given the vagaries of H a extinction cor- 
rections and the H2 - CO conversion factor dKennicut j|l998|; 

Komugi et ail 120051 ; iKennicutt et. al.l 



Wong fc Blit3 2002: 



20071 ; iBoissier et. ail 



20071 '). Also, the initial metallicity in 



model J must be non-zero, otherwise there are too many 
metal-poor stars. This model does a poor job of repro- 
ducing the observed present-day gas surface density of ~ 
0.8 M Q pc -2 (see ©. The gas density in model J decreases 
from about 3 to 2.3 Mq pc -2 for the Padova tracks and from 
about 2.2 to 1.5 Mq pc -2 for the Teramo tracks. Hence, this 
quantity is above the range plotted in panel (a) of Fig. [9] In 
model K, the AMR decreases over the last 2 Gyr because gas 
inflow continues while very few stars form. Model K also has 
more gas inflow (~ 20%) occurring in the last 3 Gyr than 
the original free inflow model (< 10%). 



6 COMPARISON TO OBSERVATIONS IN M33 

Next, we will examine in further detail the predictions of 
the original free inflow solutions, since they provide the best 
fits. The Padova solutions have a present-day gas mass of 
~ 0.4 ± 0.2 Mq pc~ 2 . From the difference between the 
Padova and Teramo solutions, we estimate a systematic er- 
ror of a factor of ~ 2 due to uncertainties in the stellar 
tracks. The projected mean HI column density in this field 
is ~ 1.1 Mq pc -2 , as measured from a mosaic constructed 
from Very Large Array and Green Bank Telescope observa- 
tions (D. Thilker, private communication). If we correct for 
a disc inclination of 56° and for a helium abundance one- 
third that of hydrogen, then this becomes ~ 0.8 Mq pc -2 . 
The true inclination of the HI disc at this location is some- 
what uncertain because of th e well-known warp at large radii 
l|Corbelli fc Schneidedll997h . A change of ±10° in the incli- 
nation would change the surface density by approximately 
± 0.2 Mq pc -2 . For comparison, the azimuthally-averaged 
HI column density at this r adius in M33 is ~ 3 Mq pc -2 
l|Corbelli fc Schneider|[l997l ; ICorbellilliool ). 

The present-day SFR provides another useful check 
on our results. Using GALEX near-UV and far-UV im- 
ages of M33, Boissier et al. (2007; private communication) 
computed the UV surface brightness in our fiel d and then 
conve rted that to a SFR using the relation in Ke nnicuttl 
(1998), whic h is a slightly modifi ed version of the one first 
presented in iMadau et al.l (|l998l ). The infrared data used 
to estimate extinction was of poor quality this far out in 
M33, so the SFR could only be constrained to the range 
~ (0.6 — 0.9) x 10 -4 Mq yr _1 , where the lower limit corre- 
sponds to zero extinction and the upper limit corresponds 
to an extinction of Afuv = 0.49. The azimuthally-averaged 
UV SFR at this radius lies at the lower limit. The chief 
sources of uncertainty in these limits are the UV flux-SFR 
calibration, which was based on theoretical isochrone and 
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Figure 9. Effects of varying the model parameters on the gas mass and chemical composition using the Padova tracks. The solid line 
with diamonds is the original free inflow model and the other lines are tests A through K (see text for details). Only the original free 
inflow model errors are shown for the sake of clarity. The panels are identical to Fig. [6] 



spectral libra ries, and the a ssump tion of a Salpeter IMF. 
According t o iMadau et~aH |l99&1), adopting a Scalo IMF 
l|Scalolll986h . which is more similar to the IMF we have 
used, the resulting SFR would be ~ 50% smaller. Therefore, 
the range above becomes ~ (0.3 — 0.5) x 10~ 4 Mq yr -1 , 
which is in good agreement with our predictions of ~ 0.4 
and 0.3 x 10 ~ 4 Mq yr -1 for the Padova and Teramo free 
inflow solutions, respectively. 

Another i mportant check com es from oxygen abun- 
dances in M33. iMagrini et alj (|2007T ) compiled an extensive 
catalog of previously measured abundances in Hll regions, 
type A-B supergiant stars, and planetary nebulae (PNe). 
This catalog is plotted in Fig. 1 101 with the predictions of the 
Padova and Teramo free inflow models for the present-day 



oxygen abundance (star and diamond, respectively), whose 
abscissa values are offset from each other for clarity. The H II 
regions and supergiants probe the present-day ISM abun- 
dance whereas the PNe probe the ISM abundance at older 
ages. The masses and lifetimes of the PNe progenitors are 
not known with great certainty. The sample of Magrini et al. 
covers only the brightest two magnitudes of the PNe lumi- 
nosity function, so it could be biased toward the high end of 
the progenitor mass range. Therefore, the progenitors prob- 
ably had MS lifetimes less than a few Gyr. 

Most of the H II region abundances in the Magrini et al. 
compilation were made from direct T e measurements, gen- 
erally considered the most reliable kind. Eight objects have 
at least two independent measurements and three of those 
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Figure 10. Oxygen abundances as measured in (a) Hn regions, 
(b) A-B supergiant stars, an d (c) PNe. The solid and d ashed 
lines in panel (a) co me from iRosolowskv fc Simo n (2008) and 
IViironen et all (120071 ), respectively. In all three panels, the star 
and diamond are, respectively, the Padova and Teramo free inflow 
model predictions for the present day oxygen abundance. Their 
abscissa values are offset from each other for clarity. 



objects have three independent measurem ents. We have 
also a dded the recent gradient derived by IViironen et al.l 
lj2007l ) (dashed line) for ~ 60 H II regions based on the 
log(Hq/[S ii]AA6717+ 6731) vs. log(Hq/[N u]A6583 ) diagnos- 
tic diagram. Recently. IRosolowskv fc Simon (|2008h reported 
oxygen abundances of 61 H II regions in the southwest region 
of M33 derived from the temperature sensitive emission line 
[Oiii] A4363 A. They found a shallow gradient which we 
show as the solid line in Fig. 1101 

Our results are more consistent with an overall shal- 
low gradient or a fl attening in the gradient at large radii 
l|Magrini et al.ll2007h than they are with a constant gradient 
over M33's entire disc. However, even if the gradient was flat 
beyond Rd p ~ 30', our results would be several tenths of a 
dex higher than expected from the data in Fig. 1101 The mag- 
nitude of this discrepancy is comparable to the random and 
systematic uncertainties in our results, so it may not be sta - 
tistically significant. Moreover. IRosolowskv fc Simonl (|2008l ) 
noted a scatter of 0.11 dex in their sample unattributable 
to measurement errors, possibly indicating chemical inho- 
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Figure 11. Comparing the [a/Fc] vs. [Fe/H] relation for M33, the 
SV, LMC, and SMC. The solid lines are the fre e inflow so lutions 
while the other lines are the chemical models in |PT95l and [PT9Sl . 

mogeneities in M33's ISM. This fact also highlights the im- 
portance of using large samples of objects at many position 
angles within the disc. More importantly, the data in Fig. 
1101 have their own systematic uncertainties. For example, 
the abundances of H II regions tend to be lower than super- 
giants and PNe by ~ 0.2 dex. This offset could indicate a 
problem in how the temperat ure and ionization structures of 
Hn regions are treated (e.g., Garcfa-Roias fc Estebanll2007l ; 
lEsteban et ail 12004 iKennicutt et al.l I2003I . and references 
therein), but larger samples need to be compared before the 
nature of this discrepancy can be fully understood. 



7 DISCUSSION 

7.1 Comparison to Models of the Solar Vicinity 
and Magellanic Clouds 

In Fig. 1111 we compare the [a/Fe] vs. [Fe/H] relation in this 
region of M33 to that of the SV, LMC, and SMC derived in 
|PT95l and |PT98l . The models for these other three systems 
have been cited often in the literature because of their sim- 
plicity and ability to reproduce observations fairly well. The 
solid lines show the free inflow solutions using the Padova 
and Teramo tracks without being averaged over each SFH 
bin. Viewing the relation in this manner facilitates compar- 
ison with the other systems and can lead to greater physi- 
cal insight, as long as we remember that the overall shape 
is most robust. This last fact can be appreciated from the 
overall similarity of the acceptable inflow solutions in Fig. HJ] 
Fig. [EH reveals that M33's [a/Fe] vs. [Fe/H] relation is 
generally similar to the other systems. Because we have used 
almost identica l chem ical evolut ion equations and identical 
stellar yields as |PT95l and |PT98l . our models show a similar 
overall behavior in response to gas flows. The discontinuities 
in the relation for each system can be traced back to the 
interplay between the SFR, IFR, and the delayed injection 
of iron into the ISM from SNe la. Therefore, the differences 
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between the relations arise primarily from differences in the 
SFH and IFH. 

In the IRA, the abundance ratio of any two elements 
is a constant determined by the ratio of their nucleosyn- 
thetic yields. Thus, the initial [a/Fe] ratio is log(t/ a /t/Fe) — 
log(X Q /X Fe )0 « 0.4. After 1.3 Gyr have elapsed, the DPA 
begins operating as the first SN la inject iron into the ISM 
and produce the first downturn, or "knee", in the relation. 
The [Fe/H] location of this knee is dependent on the SFR, 
which, in turn, depends on the IFR and outflow rate. The 
SFR was on average higher in the SV than in the other sys- 
tems, so it experienced a faster enrichment and the knee 
occurs at a higher metallicity. 

The large loop in the Padova model results from three 
successive steps: a sudden large increase in the IFR at an 
age ~ 6 Gyr which dilutes metals in the ISM (see the IFH in 
Fig. |4ji), a corresponding increase in the SFR which raises 
[a/Fe], and the subsequent drop in [a/Fe] 1.3 Gyr later due 
to SN la. The Teramo model experiences a much smaller 
loop because it does not have as large an increase in the 
IFR at 6 Gyr. 

7.2 Implications for the Formation of the MW's 
Halo 

The ACDM hierarchical picture of galaxy formation pre- 
dicts that dark matter haloes are built up from the accre- 
tion/merging of smaller su bhaloes similar t o the protogalac- 
tic fragments proposed by ISearle fc Zin n (1978). This pic- 
ture can be tested by comparing the chemical composition of 
the MW's stellar halo field population to such subhaloes or 
their possible present-day analogs, such as the MW's dSph 
satellites. As reviewed in detail bv lGeisler et al.l (|2007h and 
references therein, the MW's field halo population is chem- 
ically distinct from these other systems, indicating that the 
former was not built up from the latter. In particular, the 
present-day dSphs and other LG galaxies are characterized 
by [a/Fe] ratios lower than most of the MW's halo, especially 
at [Fe/H] > — 2. This is commonly interpreted to mean that 
the SFR in the dSphs was lower and, therefore, their chem- 
ical enrichment was too slow to reach the halo's metallic- 
ity with a high [a/Fe] before SNe la contributed significant 
amounts of iron. 

It is interesting to compare our results for M33's outer 
disc to the spectroscopic measurements in other LG sys- 
tems. Such an exercise places our results into the context 
of the LG as a whole and could have some implications for 
the formation of the MW's halo. To that end, in Fig. [12] 
we summarize the high-resolution spectroscopic measure- 
ments of [a/Fe] a nd [Fe/H] in these systems. This figure 
is reproduced from lGeisler et all l|2007h except that we have 
added the Padova and Teramo free inflow solutions as the 
black and magenta lines Q. The green line represents the 
halo's dissipative collapse component while the green open 
st ars represent the halo 's accretion component, as defined 
bv lGratton et all l|2003T ) based on orbital kinematics. These 
authors postulated the dissipative component formed during 
the collapse of the MW's gaseous halo, as earlier outlined by 

2 To be consistent with iGeisler et al | d2007h . we plot [a/Fe] = 
([0/Fe] + [Mg/Fe] + [Si/Fe] + [Ca/Fe] + [Ti/Fe])/5 in Fig.[TJ 



lEggen et all l| 19621) , and the accretion component originated 
in the accreted subhaloes. 

The green open diamonds at [Fe/H] ~ —3 represent 
extrem ely metal-poor halo giants analysed bv lCavrel et al.l 
for which no separation into kinematic components 
has been made. All other green open points represent sam- 
ples of MW halo stars selected for their unusual orbital 
properties and which therefore are additional candidates for 
originating in chaotically accreted subhaloes. The other sys- 
tems represented in Fig. [12] are low-mass dSph satellites of 
the MW (blue filled circles, stars, and asterisks), Sagittarius 
dSph (red filled circles and stars), LMC (cyan filled squares 
and upward pointing triangles), and several dlrrs (orange 
filled downward pointing triangles). Typical random mea- 
surement uncertainties are i ndicated in the bott om left cor- 
ner. We refer the reader to IGeisler et alj (|2007T ) for a more 
detailed explanation and for a list of all the references for 
the data. 

Based on our results, if we could obtain high resolu- 
tion spectroscopic abundances for the location we have stud- 
ied in M33, we would expect most of its stars to lie near 
([a/Fe]) ~ 0.1 and ([Fe/H]) ~ -0.7. Of all the systems pre- 
sented in Fig. 1121 this region of parameter space is most 
similar to the LMC. We also note that M33 resembles the 
MW's halo more than the low-mass dSphs. The plateau at 
[Fe/H] < —2.0 was built into the models via the chemical 
yields (see t|7.1|l . but the key fact is that M33 maintains a 
higher [a/Fe] ratio at a given [Fe/H] relative to the low- 
mass dSphs. The agreement between M33 and the MW's 
halo becomes progressively worse at [Fe/H] > —2. However, 
M33's relation overlaps with some of the candidate accre- 
tion stars suggesting they could have originated in a system 
whose SFH, AMR, and IFH were similar to M33's outer 
disc, though not necessarily a spir al galaxy like M33. T his is 
consistent with the suggestion of IGeisler et al.l (|2007) that 
some of the accreted MW halo stars originated in high mass 
dwarf systems. 

Although we have little information on M33's in- 
ner disc evolu t ion, t he presence of a metallicity gradient 
( Mag rTni et alJ | 2007h and a wide v ariety of ages inferred 
from CMDs (|Saraiedini et all 120001 ) hint at an extended 
SFH with a faster enrichment taking place there. There- 
fore, we hypothesize that the knee in the [a/Fe] vs. [Fe/H] 
relation in M33's inner disc occurs at a higher metallic- 
ity than in the outer disc. This brings us to another key 
point, that candidate accretion halo stars with different 
[a/Fe] but the same [Fe/H] did not necessarily originate in 
different objects. In other words, the a-element and iron 
abundances alone do not necessarily provide enough infor- 
mation to tag individual halo star s to unique progenitors 
l|Freeman fc Bland-Hawthorrj|20o3) ■ This is because the in- 
ner regions of a protogalactic fragment could have become 
enriched with metals faster than the outer regions. There- 
fore, the candidate accretion stars that we observe today 
with low and high [a/Fe] could have originated in the outer 
and inner regions of the fragment, respectively. 

How realistic is such a scenario? How much time would 
have been necessary and available for this hypothetical frag- 
ment to attain the range of [Fe/H] and [a/Fe] observed today 
in the MW halo's accreted stars? Our results indicate M33's 
outski rts took ~ 4 Gyr to reach [Fe/H] > —1.0 and [a/Fe] 
~ O.l. lZentner fc Bullock! l|2003h used a semi-analytic model 
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Figure 12. Comparing the [o/Fe] vs. [Fe/H] relation in one location in M33's outskirts to spectro scopic a bunda nce ratios of individual 
stars in various LG systems. The data were compiled from the literature and presented in lGeisler et al.l [I2007TI . to which we refer the 
reader for more details. The systems shown are the MW halo (green open points and green line), low-mass dSphs (blue filled circles, 
stars, and asterisks), Sagittarius dSph (red filled circles and stars), LMC (cyan filled squares and upward pointing triangles) and dlrrs 
(orange filled downward pointing triangles). The Padova and Teramo free inflow models are black and magenta lines, respectively. The 
error bars in the lower left corner represent typical random measurement uncertainties. 



to examine the survival of orbiting subhaloes in the MW's 
potential. They found that for a typical subhalo orbit, the 
survival time after entering the MW's halo ranged from 3 
to over 8 Gyr, depending on the subhalo's mass and con- 
centration relative to the MW's. Any time spent outside the 
MW halo prior to accretion would only relax the necessary 
post-accretion survival time. Thus, it seems entirely plausi- 
ble that this hypothetical fragment could have had enough 
time to develop a population gradient before disruption. 



Present-day popul ation gradients have been found in 



most of the dSphs (e.g., Harbcck et al.ll 


2001 ; Winnickl 


2003]; 


Tolstoy et all 12004 iKodi et al.1 l200fj; 
Bellazzini et al.l 120051; Battaelia et al. 


Stetson et al.l 
120061; fFaria 


1998|; 
et al. 


20071; iKomivama et al.ll2007l; McConnachie et al.l 20071). in 



120071). NGC 6822 (Ide Blok fc Walter} I2006I ). and in M33 

In most cases, the 



( Rowe et al] 120051. [Paper ij , iPaper 
inner regions of the galaxies are more metal rich and/or 
older than the out er regions. In another recent study, 
iBernard et al.1 (2008) found that RR Lyrae stars in the cen- 
tral region of Tucana were more metal-rich than in the outer 
region, suggesting a metallicity gradient was already present 
in the first few Gyr of that galaxy's history. 

The evolution of metallicity gradients in disc galaxies 
is a matter of some debate, with some studies predicting 
the gradient to increase with tim e and others predicting 
the opposite (|Magrini et al.l [20071 . and refer ences therein). 
The evolutionary model of M33 presented in iMagrini et al.l 



l|2007l ) predicts the metallicity gradient to have been steeper 
in the past. Measurements of the spatially resolved SFH and 
AMR of these low mass galaxies could help us understand if 
they can develop chemical gradients in the first few Gyr of 
their lives or if such gradients were present in the gas from 
which they formed. 



8 CONCLUSIONS 

In the framework of a simple chemical evolution scenario 
which adopts instantaneous and delayed recycling for the 
nucleosynthetic products of Type II and la SNe, we have 
modelled the observed CMD in one location in M33's out- 
skirts. In this scenario, interstellar gas forms stars at a rate 
modulated by the KS relation and gas outflow occurs at a 
rate proportional to the SFR. Compared to the common 
CMD fitting method of allowing age and metallicity to be 
free parameters, this scenario yields a more physically self- 
consistent SFH with fewer free parameters and makes more 
predictions that can be tested with independent observa- 
tions. Moreover, this method puts broad constraints on the 
role of gas flows and extends the method of chemical finger- 
printing to stellar systems, like M33, that are beyond the 
reach of current high resolution spectrographs. 

The precise details of our results depend on which stel- 
lar evolutionary tracks are used, Padova or Teramo. Nev- 
ertheless, the broad trends appear to be robust. We found 
that, when the star formation efficiency, e, is constant, the 
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canonical closed box model fails to reproduce the observed 
distribution of stars in the CMD. Models with an expo- 
nentially declining IFR from 14.1 Gyr ago exhibit similar 
discrepancies, but to a lesser magnitude. Instead, the in- 
flow models which best reproduce the observed CMD have 
a significant fraction (> 50%) of gas inflow taking place in 
the last 7 Gyr and a smaller fraction (< 10%) taking place 
within the last 3 Gyr. This lea ds to a SF H in overall agree- 
ment with what we found in IPaper III and suggests the 
traditional method of synthetic CMD fitting can be physi- 
cally self-consistent. 

Allowing e to vary with time, as might be expected if 
the star formation time-scale depends on local ISM proper- 
ties, significantly improves the closed box model if the initial 
metallicity is [Fe/H] ~ —1.3, but the resulting present-day 
gas mass surface density is too high compared to the ob- 
served value. A model with a varying e and constant, non- 
zero IFR reproduces the CMD and gas mass about as well 
as the best constant e models. In this case, more inflow can 
take place in the last 3 Gyr, but there is still a significant 
fraction (~ 50%) of inflow occurring over the last 7 Gyr. The 
amount of variation in e required by these models, however, 
could be larger than the intrinsic scatter in the KS relation 
of other galaxies. 

We also examined the [a/Fe] vs. [Fe/H] relation, a key 
diagnostic of the evolutionary history of stellar systems. Like 
the MW's dwarf satellites, the bulk of stars at this location 
in M33's outskirts have [a/Fe] ratios lower than the MW's 
halo field. Stars formed over 8 Gyr ago have ([a/Fe]) ~ 0.2± 
0.1 and stars forming today have {[a/Fe]} ~ —0.1 ±0.1. The 
mean [a/Fe] ratio of all stars ever formed at this location in 
M33 was found to be ~ 0.1 ± 0.1. From the tests we have 
conducted, we estimate that the systematic errors of these 
values are ~ ±0.2 dex, comparable to that of some high 
resolution spectroscopic measurements of other systems (e.g. 
iMonaco et aT]|2005l ; iGeisler et al1l2007l ). 

Our results paint a picture in which M33's outer disc 
formed from the protracted inflow of gas over several Gyr 
with at least half of the total inflow occurring since z ~ 1, 
but relatively little since z ~ 0.25. All the acceptable in- 
flow models have a similar IFR ~2x 10 -4 Mq yr _1 kpc -2 
averaged over the lifetime of the Universe. This estimate is 
uncertain by at least a factor of 2 and the IFR could have 
intrinsic variations of a factor of ~ 3 when averaged over 2 
— 3 Gyr time-scales. Nevertheless, it still gives a rough indi- 
cation of the mean gas IFR that could occur in the outskirts 
of low mass spiral galaxies. 

A useful baseline for co mparison is the li fetim e average 
IFR in the SV, predicted bv lColavitti et al.l (|2008l ) to be <~ 
4 x 10~ 3 Mq yr _1 kpc -2 . This estimate comes from several 
chemical evolution models which adopt different shapes for 
the IFH and which reproduce many local observables, like 
the abundances of various elements, MDF of long-lived stars, 
and present-day gas fraction, total mass density, SFR, IFR, 
and SN rate. That the mean IFR in the SV has been higher 
than in M33's outer disc is not surprising since the MW is 
more massive than M33 and the SV is at a much smaller 
radius (in terms of disc scale lengths) than the M33 field we 
hav e studied. 

ISommer-Larsen et all l|2003h ran an N-body cosmolog- 
ical simulation and selected 12 dark matter haloes for more 
detailed simulations which included baryons. Two Milky 



Way-type spirals, SI and S2, that resulted from the sim- 
ulations experienced gas accretion at a mean rate of ~ 
10~ 3 Mq yr -1 kpc -2 at radii of 6 — 7 disc scale lengths. 
This rate is about 5 times larger than the mean IFR we 
have found at a similar location in M33's outer disc, but 
SI and S2 were about 5 times more massive than M33. 
The total disc averaged accretion rates of the least mas- 
sive simulated galaxies, with rotation velocities comparable 
to M33's, were ~ 15 — 35% smaller than those of SI and S2 
at z = 0. Therefore, if the outer disc accretion rate scales 
with a galaxy's rotation velocity in the same way as the 
total disc averaged rate, and if this scaling holds for all red- 
sh ifts, then our results are con sistent with the simulations 
of ISommer-Larsen et al.l (|20031 ) . 

What is the source and nature of the gas inflow? If 
it occurs in the disc plane, it could be caused by spi- 
ral density waves, viscosity in the disc gas, or gas with 
lower angular momentum falling onto the disc at larger 
radii and flowing toward M33's nucleus ( Lacev fc Falllll985l; 



Lin fc Pringielll987 
200d ; iRoskar et al.l 



.Bertin fc Linlll996l ; IPortinari fc Chiosil 
20081 ). Gas flowing from above the disc 



plane could come from the condensation of a hot halo corona 
as described in |T] However, such a process is expected to 
be more efficient in a ma ssive spiral like the MW than in a 
late- type spiral like M33 (jDekel fc Birnboiml [20061 '). 

The condensed clouds predicted by numerical simu- 
lations to fall onto a disc galaxy have properties simi- 
lar to the MW's HVC population (|Peek et all 120081 ) . Re- 
cent surveys of HI emission around M31 and M33 have re- 
vea led an analogous HVC population around these galax- 
ies dThilker et al.ll2004l; iBraun fc Thilker|[2(5oi ; iGrossi etHI 
120081 ). iThilker et all (|2004l ) detected 20 clouds within 50 kpc 
of M31's disc and the maps presented by I Gross! et al.l ( 2008) 
show many within ~ 20 kpc of M33 with an inferred total 
mass > 5 x 10 7 M Q . The M31/33 clouds appear to have a 
mean mass ~ 10 5 M^ and size ~ 1 kpc (IWestmeier et all 
l2005l ; lGrossi et al.1120081 ). The lifetime integrated inflow mass 
of ~ 10 6 Mq in the M33 region we have studied requires 
the accretion of ~ 10 such clouds. Extrapolation to M33's 
entire disc is risky given how little we know about its evolu- 
tion and the population of clouds around M33. In addition 
to condensed halo gas, the clouds could also be low-mass 
dark matter subhaloes that never formed sta rs or the tidal 
debris from recent mer gers and interactions (jThilker et al.l 
2004; Grossi et al.ll2008T ). In the future, we plan to apply the 
techniques presented in this paper to other regions of M33 
and other galaxies. This approach can provide insights on 
the origin and amount of accreted gas and coupling (or lack 
thereof) between the baryonic and dark matter accretion 
histories of these systems. 
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